**********************
*** THIS SCRIPT PROVIDES SUPPORTING CALCULATIONS FOR NUMBERS PROVIDED IN THE ONLINE APPENDIX AND MAIN TEXT
*** IT ALSO CREATES THREE APPENDIX TABLES: TWO VALIDATING ADMINISTRATIVE PRODUCTIVITY VARIABLES, and ONE DIFFERENTIAL ATTRITION TABLE
**********************
clear 
set more off
local texsave_settings "replace autonumber nofix"

* All required add-ons are stored in the /packages and /auxiliary folders
adopath ++ "$Wellness_WhatDoesWWDo/scripts/packages"
adopath ++ "$Wellness_WhatDoesWWDo/scripts/auxiliary"
mata: mata mlib index

**********************************************
* Table reporting the loadings on the productivity indices (principle component analysis)
**********************************************

* Loadings
tempfile yr0 yr1 yr2

foreach v in yr0 yr1 yr2 {
	use "$Wellness_WhatDoesWWDo/results/intermediate_files/productivity/prod_index_`v'_loadings.dta", clear
	cleanvars var
	ren loadings_`v'_1 loadings_`v'
	save "``v''", replace
}

use `yr0', clear
merge 1:1 var using "`yr1'", nogenerate
merge 1:1 var using "`yr2'", nogenerate

gen nonmissing = mi(loadings_yr0)
order var
sort nonmissing loadings_yr0 loadings_yr1
drop nonmissing

tostring loadings_*, format(%5.3fc) force replace
replace loadings_yr0 = "N/A" if loadings_yr0=="."
replace loadings_yr1 = "N/A" if loadings_yr1=="."
replace loadings_yr2 = "N/A" if loadings_yr2=="."

label var var "Productivity Variables"
label var loadings_yr0 "Baseline"
label var loadings_yr1 "First-Year Follow-Up"
label var loadings_yr2 "Longer-Run Follow-Up"
local title "The Loadings of the First Principal Component of Productivity"

local fn `"Notes: The first principal component of productivity corresponds to the LEFTQUOTES_productivity index'' outcome reported in other tables. This component is calculated as a linear combination of productivity variables, where the weights in that linear combination are equal to the loadings reported in this table and the variables in that linear combination are normalized to have a mean of zero and variance of one. The sum of the squared loadings is equal to one for each column."'
texsave using "`yr0'", `texsave_settings' title("`title'") varlabels footnote(`"`fn'"', size(footnotesize)) marker(tab:appendix_productivity_loadings)
filefilter "`yr0'" "$Wellness_WhatDoesWWDo/results/tables/appendix_productivity_loadings.tex", from(`"LEFTQUOTES_"') to("`=char(96)'`=char(96)'") replace


**********************************************
* Productivity validation table
**********************************************

tempfile validation_0717 tmp

local job_vars_admin_0717  "salaryRaise_0616_0717 promotion_0616_0717 titleChange_0616_0717 terminated_0717 sickleave_0816_0717"
local job_vars_survey_0717 "sickdays_0717 hrsworked50_0717 jobsatisf1_0717 jobsatisf2_0717 mgmtsafety_0717 happywork_0717 productive_0717 promotion_0717 jobsearch1_0717 jobsearch2_0717"

use "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", clear

* "The average salary raise in our sample was 5.9 percent after one year"
summ salaryRaise_0616_0717 if !mi(treat)
assert abs(r(mean)-0.059)<0.0005

* "For those with a job title change, the average raise was 14.5 percent."
summ salaryRaise_0616_0717 if !mi(treat) & titleChange_0616_0717
assert abs(r(mean)-0.145)<0.0005

* "A small number ($<5$ percent) of employees with job title changes did not receive an accompanying salary raise"
count if titleChange_0616_0717==1 & !mi(treat) & !mi(salaryRaise_0616_0717)
assert r(N)==763
count if titleChange_0616_0717==1 & !mi(treat) &     salaryRaise_0616_0717 <0.0001
assert r(N)==37
assert 37/763<0.05

* Calculate (constrained) outcome means
gen to_use = 1
foreach v of varlist `job_vars_survey_0717' {
	replace to_use = 0 if mi(`v')
}

foreach v of varlist `job_vars_admin_0717' {
	summ `v' if to_use==1
	local mean_`v' = string(r(mean),"%5.3f")
	local N_`v' = string(r(N),"%5.0fc")
}

preserve
local run_no = 0
qui foreach xvar of varlist `job_vars_survey_0717' {

	local replace replace
	foreach y of varlist `job_vars_admin_0717' {
	
		reg `y' i.`xvar', robust
		regsave 1.`xvar' using `tmp', t p table("`y'", asterisk(10 5 1) format("%5.3f")) addlabel(outcome,"`y'", regressor,"`xvar'") `replace'
		local replace append
	}
	
	use "`tmp'", clear

	gen run_no = `run_no'
	gen outcome = "`xvar'"
	keep if strpos(var,"_coef") | strpos(var,"_stderr") | strpos(var,"_pval") | inlist(var,"N")
	replace var = "zzN" if var=="N"
	replace var = subinstr(var,"1.`xvar'","`xvar'",.)	
	
	if `run_no'>0 append using "`validation_0717'"
	save "`validation_0717'", replace		
	local run_no = `run_no'+1
	restore, preserve
}
restore, not

use "`validation_0717'", clear
sort run_no var

set obs `=_N+2'
replace var = "Sample size (outcome mean)" in `=_N-1'
replace var = "Outcome mean" in `=_N'

foreach v of varlist *_0717* {
	replace `v' = "("+`v'+")" if strpos(var,"_stderr")
	replace `v' = "\textit{N=}" + `v' if var=="zzN"
	
	replace `v' = "`mean_`v''" if var=="Outcome mean"
	replace `v' = "`N_`v''" if var=="Sample size (outcome mean)"
}
replace var = "" if strpos(var,"stderr")
replace var = subinstr(var,"_coef","",.)
drop if strpos(var,"_pval")
cleanvars var
ingap 3 6 9 12 15 18 21 24 27 if var=="zzN", after
replace var = "" if var=="zzN"
drop run_no outcome

label var salaryRaise_0616_0717 "Annual salary (share of baseline salary)"
label var promotion_0616_0717 "Job promotion"
label var titleChange_0616_0717 "Job title change"
label var terminated_0717 "Job terminated"
label var sickleave_0816_0717 "Sick leave (days/year)"

local title "Associations Between Administrative and Survey Measures of Productivity"
local fn `"Notes: Each row and column reports estimates from a separate regression, where observations include all individuals who completed the one-year follow-up survey. The outcome in each regression is specified by the table column. The independent variable is the survey response, which is always an indicator variable. (The Presenteeism survey variable is omitted from this table.) The outcome mean is calculated for the sample of observations that have non-missing values for all survey variables listed in this table. Regressions are unweighted. Robust standard errors are reported in parentheses. A */**/*** indicates significance at the 10/5/1\% level using conventional inference, i.e., not adjusting for multiple outcomes."'
texsave using "$Wellness_WhatDoesWWDo/results/tables/appendix_productivity_validation.tex", `texsave_settings' title("`title'") hlines(-2) size(scriptsize) varlabels footnote(`"`fn'"', size(scriptsize)) marker(tab:appendix_productivity_validation)

**********************************************
* Differential attrition table
**********************************************
use "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", clear
tempfile attrition_tbl

keep if !mi(treat)
preserve

* "the difference is statistically significant at the two-year mark ($p=0.036$)"
reg survey2018_c treat, robust
matrix tmp = r(table)
matrix a = tmp["pvalue",1]
scalar p = a[1,1]
assert abs(p-0.036)<.0005

count if treat==1
local sample_size_treat = string(`r(N)',"%12.0fc")
count if treat==0
local sample_size_control = string(`r(N)',"%12.0fc")

local run_no = 0
qui foreach v in survey2018_c survey2017_c terminated_0119 terminated_0717 covg_0816_0119 covg_0816_0717 {
	
	* Calculate means for treat and control
	reg `v' ibn.treat, nocons robust
	
	* Test for equality of means
	test i0.treat==i1.treat
	local pval = `r(p)'
	
	regsave
	keep var coef
	
	gen outcome="`v'"
	replace var = "Control" if var=="0bn.treat"
	replace var = "Treatment" if var=="1.treat"
	reshape wide coef, i(outcome) j(var) string
	gen pval = `pval'
	
	if `run_no'==1 append using `attrition_tbl'
	save `attrition_tbl', replace
	local run_no = 1
	restore, preserve
}

restore, not
use `attrition_tbl', clear

ren coefControl Control
ren coefTreatment Treatment
gen Difference = Treatment - Control
format Control Treatment Difference pval %5.3f
order outcome Control Treatment Difference pval
tostring Control Treatment, format(%12.3fc) replace force

set obs `=_N+1'
replace outcome = "Sample size" if mi(outcome)
replace Control = "`sample_size_control'" if outcome == "Sample size"
replace Treatment = "`sample_size_treat'" if outcome == "Sample size"

replace outcome = "2017 survey" if outcome == "survey2017_c"
replace outcome = "2018 survey" if outcome == "survey2018_c"
replace outcome = "Health insurance enrollment (first 12 months)" if outcome=="covg_0816_0717"
replace outcome = "Health insurance enrollment (first 30 months)" if outcome=="covg_0816_0119"
replace outcome = "Job terminated (first 12 months)" if outcome=="terminated_0717"
replace outcome = "Job terminated (first 30 months)" if outcome=="terminated_0119"

label var Control "Control"
label var Treatment "Treatment"
label var Difference "Difference"
label var pval "\(p\)-value"

local fn `"Notes: This table reports health insurance enrollment rates and job termination rates for the first 12 and the first 30 months following randomization, and completion rates for the 2017 and 2018 online surveys. An individual's insurance enrollment is defined as the number of covered months divided by the length of the sample period (either 12 or 30 months). Columns (1)-(2) report unweighted means for the control and treatment groups. Column (3) reports the difference between the two means and column (4) reports the \(p\)-value from a joint test of equality of the two means."'
texsave using "$Wellness_WhatDoesWWDo/results/tables/appendix_attrition.tex", `texsave_settings' hlines(-1) title("Tests of Differential Attrition Between Control and Treatment Groups") varlabels footnote(`"`fn'"', size(footnotesize)) marker(tab:appendix_attrition)


***********************
* Reported p-values
***********************

* Selection analysis
use "$Wellness_WhatDoesWWDo/results/intermediate_files/selection_results.dta", clear

* "This pattern of advantageous selection is strongly significant using conventional inference (p=0.027)"
assert abs(pval_hra_c-0.027)<0.001 if outcome=="spend_0715_0716"

* "and remains marginally significant even after adjusting for the five outcomes in this family (family-wise $p=0.082$)"
assert hra_c=="0.082" if var=="spend_0715_0716_wypval"

* "at baseline, average annual medical spending among participants was \$1,384 less than among non-participants"
* "employees who completed the screening and HRA spent, on average, \$115.3 per month less on health care in the 13 months prior to our study than employees who did not participate"
assert hra_c=="-115.3**" if var=="spend_0715_0716_coef" & (115.3*12-1384)<0.5

* "employees participating in our wellness program were \textit{more} likely to have non-zero medical spending at baseline than non-participants, by about 5 percentage points (family-wise $p\leq 0.02$)"
assert activity_f_c=="0.049***" if var=="nonzero_spend_0715_0716_coef"
assert activity_s_c=="0.046***" if var=="nonzero_spend_0715_0716_coef"
assert real(activity_f_c)<=0.020001 if var=="nonzero_spend_0715_0716_wypval"
assert real(activity_s_c)<=0.020001 if var=="nonzero_spend_0715_0716_wypval"

* "participants in the treatment groups with \$100 and \$200 screening incentives spent, on average, \$79 \textit{more} per month ($p=0.06$)"
use "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", clear
gen group_bc = inlist(StudyArm_ScreenHRA,2,3)
regress spend_0715_0716 group_bc if treat == 1 & hra_c == 1 [aweight=covg_0715_0716], robust
test group_bc
assert abs(r(p) - 0.06) < 0.001

* "as we increase activity incentives, the marginal participant has significantly \textit{lower} spending ($p=0.03$)"
gen group_75 = StudyArm_Activities == 2
regress spend_0715_0716 group_75 if treat == 1 & hra_c == 1 [aweight=covg_0715_0716], robust
test group_75
assert abs(r(p) - 0.03) < 0.001

* Balance table: "This is confirmed by our joint balance tests, which fail to reject the null hypothesis that the variables in Panel B ($p = 0.821$), Panel C ($p = 0.764$), or Panel D ($p = 0.752$) are not predictive of group assignment"
insheet using "$Wellness_WhatDoesWWDo/results/tables/balance_tests1.tex", delimiter(&) clear
assert v5=="0.821" if v1=="Joint balance test for panel B (\(p\)-value)"

insheet using "$Wellness_WhatDoesWWDo/results/tables/balance_tests2.tex", delimiter(&) clear
assert v5=="0.764" if v1=="Joint balance test for panel C (\(p\)-value)"
assert v5=="0.752" if v1=="Joint balance test for panel D (\(p\)-value)"

* "the difference in average spending between treatment and control was only \$10.8 per month"
use "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", clear
reg spend_0816_0717 treat [aw=covg_0816_0717 ], robust
assert _b[treat]-10.8<0.05


***********************
* Reported confidence-interval values
***********************
local y terminated_0119
estimates use "$Wellness_WhatDoesWWDo/results/intermediate_files/0816_0119/treat_effects/estimates/ITT_Mean_`y'.ster"
local y_mean = _b[rhs]
estimates use "$Wellness_WhatDoesWWDo/results/intermediate_files/0816_0119/treat_effects/estimates/ITT_Lasso_`y'.ster"
regsave, ci
keep if var == "rhs"
assert _N==1
local y_ci_upper = ci_upper[1]

* "We can rule out a positive retention effect of 2.4 percentage points (12.0 percent) for iThrive."
assert abs(round(`y_ci_upper'*100, 0.1) - 2.4) < 0.000000001
assert abs(round(`y_ci_upper'*100/`y_mean', 0.1) - 12.0) < 0.000000001

* "we show that our OLS (non-RCT) estimate for medical spending... is ruled out by the 99 percent confidence interval of our IV (RCT) estimate."
local y spend_0816_0717
estimates use "$Wellness_WhatDoesWWDo/results/intermediate_files/0816_0717/treat_effects/estimates/IV_ztreat_Lasso_`y'.ster"
regsave, ci level(99)
keep if var == "rhs"
assert _N==1
local IV_ci99_lower = ci_lower[1]
estimates use "$Wellness_WhatDoesWWDo/results/intermediate_files/0816_0717/treat_effects/estimates/OLS_Lasso_`y'.ster"
regsave, ci
keep if var == "rhs"
assert _N==1
local OLS_coef = coef[1]
di "OLS coef: " `OLS_coef'
di "IV lower bound 95% CI: "`IV_ci99_lower'
assert `OLS_coef' < `IV_ci99_lower'

***********************
* Costs of the wellness program
***********************
use "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", clear

* Participation and cost analysis only done for those in the treatment group
keep if treat==1

***
* Average variable costs for screening, hra, and fall/spring wellness activities
***
* Calculate variable costs for each person in the treatment group
*  Average screening costs are $67.34
*  Average HRA costs are $10.88
*  Average wellness course costs are $26.07

* "The total variable cost for a health screening and HRA is $78.22"
local screen_cost 67.34
local hra 10.88
local activity 26.07
assert abs(`screen_cost'+`hra' - 78.22)<0.01

gen reward_hra = 100*(StudyArm_ScreenHRA=="B":StudyArm_ScreenHRA) + 200*(StudyArm_ScreenHRA=="C":StudyArm_ScreenHRA)
gen reward_activity = 25*(StudyArm_Activities=="25":StudyArm_Activities) + 75*(StudyArm_Activities=="75":StudyArm_Activities)

* Rewards only paid to those who complete. Activity costs incurred for anybody who registers, even if they don't go on to complete.
gen varcost_screen = `screen_cost'*screening_c + (`hra' + reward_hra)*hra_c
gen varcost_activity_f = `activity'*activity_f_r + reward_activity*activity_f_c 
gen varcost_activity_s = `activity'*activity_s_r + reward_activity*activity_s_c

gen varcost = varcost_screen + varcost_activity_f + varcost_activity_s

* Average variable cost among paticipants, including expected reward payments, is $271
summ varcost if hra_c==1 
assert abs(r(mean)-270.63) < 0.01


***********************
* College Wellness Programs
***********************


* 370/548=68% of colleges have a wellness program
use "$Wellness_WhatDoesWWDo/data/raw/college_wellness_programs/college_wellness_programs.dta", clear

assert _N==548
count if wellness=="Yes":yesno
assert r(N) == 370

***********************
* Baseline survey stats
***********************

* Baseline survey invitations sent to 12,459 people
use "$Wellness_WhatDoesWWDo/data/raw/HR_Baseline.dta", clear
assert _N==12459

* 49 distinct colleges, 323 distinct organizations
bysort HomeCollege: gen count1 = _n==1
count if count1

bysort HomeOrg: gen count2 = _n==1
count if count2

* 7,468 people clicked on survey, 4,918 began the survey, and 4,834 completed it (note: excludes withdrawals)
use "$Wellness_WhatDoesWWDo/data/raw/BaselineSurvey-2016.dta", clear
count if Clicked_BaselineSurvey ==1
assert r(N)==7468

count if Completed_BaselineSurvey ==1
assert r(N)==4834 

count if FinalDisposition_x10000=="Complete":FinalDisposition_x10000 | FinalDisposition_x10000=="Partial":FinalDisposition_x10000 | FinalDisposition_x10000=="Logged on to survey, not past consent language":FinalDisposition_x10000
assert r(N)==7468 

count if FinalDisposition_x10000=="Complete":FinalDisposition_x10000 | FinalDisposition_x10000=="Partial":FinalDisposition_x10000
assert r(N)==4918

count if FinalDisposition_x10000=="Complete":FinalDisposition_x10000
assert r(N)==4834 

* Out of 4,380 people who answered the age question, 24 were off by more than one year when compared to the HR data
keep if Completed_BaselineSurvey 
keep AnalysisID Background_Gender Background_Age Time*
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/raw/HR_Baseline.dta", keep(match) assert(match using) nogenerate keepusing(Gender BirthYear)

assert inlist(Gender,"M","F")
tab Background_Gender if Gender=="M"
tab Background_Gender if Gender=="F"

count if !mi(Background_Age)
assert r(N)==4830

gen age_HR = 2016 - BirthYear
assert !mi(age_HR)
gen diff = age_HR - Background_Age

count if abs(diff)>1 & !mi(Background_Age)
assert r(N)==24

* Quickest completion time: 3.2 minutes
gen time_min = (TimeCompleted - TimeFirstClick )/(1000*60)
summ time_min

* Conditional on finishing in less than one hour, average completion time was 15 minutes
summ time_min if time_min<60
assert abs(15-r(mean))<0.1

* Out of the pre-registered set of 21 variables, 13 people are missing at least one response
local svy_hvars "everscreen active active_try cursmk othersmk formsmk drink drinkhvy chronic health1 health2 problems energy ehealth overweight badhealth sedentary"
local svy_uvars "druguse physician hospital"
local productivity_vars "sickdays"
use "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", clear
drop if StudyArm=="Not Randomized":StudyArm
gen miss=0
gen count = 0
foreach v in `svy_hvars' `svy_uvars' sickdays {
	replace miss=1 if mi(`v')
	replace count = count+1
}
count if miss==1
assert r(N)==13
assert count==21
drop miss count

**************************
* Size of treatment groups
**************************
* 1,534 assigned to control
use "$Wellness_WhatDoesWWDo/data/raw/TreatmentAssignment-2016.dta", clear
count if StudyArm=="Control":StudyArm
assert r(N)==1534

* Treatment group sizes range from 548-552
collapse (count) sum=AnalysisID if Randomized==1 & StudyArm!="Control":StudyArm, by(StudyArm) fast
sort sum
assert sum==548 if _n==1
assert sum==552 if _n==_N

assert sum==551 if StudyArm=="A25":StudyArm
assert sum==549 if StudyArm=="A75":StudyArm
assert sum==552 if StudyArm=="B25":StudyArm
assert sum==548 if StudyArm=="B75":StudyArm
assert sum==551 if StudyArm=="C25":StudyArm
assert sum==549 if StudyArm=="C75":StudyArm

**************************
* Size of the strata
**************************
use "$Wellness_WhatDoesWWDo/data/raw/TreatmentAssignment-2016.dta", clear
keep if Randomized==1 

* 69 strata
bysort Strata: gen s_count = _n==1
count if s_count==1
assert r(N)==69

* Sample size in each strata ranges from 20 to 251
gen dummy=1
collapse (sum) dummy, by(Strata) fast
summ dummy
assert r(min)==20
assert r(max)==251

**************************
* May 30, 2017 employment data
**************************
use "$Wellness_WhatDoesWWDo/data/raw/Employment-2017-08-16.dta", clear
assert _N==12459
count if TerminationDate<=mdy(8,15,2017)
assert r(N)==1537


**************************
* Health screening data
**************************

* 2016: 2047 signups 
use "$Wellness_WhatDoesWWDo/data/raw/AcuityScheduling-2016.dta", clear
count if !mi(DateScheduledFirst)
assert r(N)==2047

* 2016: 1900 screened
use "$Wellness_WhatDoesWWDo/data/raw/Biometrics_DigitizedByRAs-2016.dta", clear
count if CompletedScreening==1
assert r(N)==1900

* 2017: "In total, 2,409 employees were assigned to the \$0 follow-up screening incentive, while 2,410 employees were assigned to the \$125 follow-up screening incentive."
use "$Wellness_WhatDoesWWDo/data/raw/TreatmentAssignment-2017.dta", clear
assert !mi(StudyArm_2017)
drop if StudyArm_2017 =="Not Randomized":StudyArm_2017
count if inlist(StudyArm_2017,"T0":StudyArm_2017,"X0":StudyArm_2017)
assert r(N)==2409
count if inlist(StudyArm_2017,"T125":StudyArm_2017,"X125":StudyArm_2017)
assert r(N)==2410

* 2017: "A total of 2,004 people from both groups were successfully screened"
use "$Wellness_WhatDoesWWDo/data/raw/Biometrics_DigitizedByRAs-2017.dta", clear
count if CompletedScreening==1
assert r(N)==2004

* 2018: "In total, 2,399 employees were assigned to the \$0 follow-up screening incentive, while 2,400 employees were assigned to the \$75 follow-up screening incentive."
use "$Wellness_WhatDoesWWDo/data/raw/TreatmentAssignment-2018.dta", clear
assert !mi(StudyArm_2018)
drop if StudyArm_2018 =="Not Randomized":StudyArm_2018
count if inlist(StudyArm_2018,"T0":StudyArm_2018,"X0":StudyArm_2018)
assert r(N)==2399
count if inlist(StudyArm_2018,"T75":StudyArm_2018,"X75":StudyArm_2018)
assert r(N)==2400

* 2018: "A total of 2,004 people from both groups were successfully screened"
use "$Wellness_WhatDoesWWDo/data/raw/Biometrics_DigitizedByRAs-2018.dta", clear
count if CompletedScreening==1
assert r(N)==1761


**************************
* Online health assessment data
**************************

* 1,848 people completed the HRA
use "$Wellness_WhatDoesWWDo/data/raw/HRA-2016.dta", clear
count if CompletedHRA
assert r(N)==1848

* Conditional on finishing in less than one hour, average completion time was 14 minutes
gen time_min = (FirstCalcDate - AssmtStarted)/(1000*60)

summ time_min if time_min<60
assert abs(13.92-r(mean))<0.1

**************************
* Wellness Activities data
**************************

* 903 completed the Fall 2016 wellness activities (out of 1304)
* "Nearly 80 percent of people who registered were signed up for HealthTrails"
use "$Wellness_WhatDoesWWDo/data/raw/Activities-2018.dta", clear
count if CompletedActivityFall16
assert r(N)==903

count if RegisteredActivityFall16
assert r(N)==1304

count if ActivityF16=="HealthTrails"
assert (r(N)/1304 -.788)<0.01

* "In total, 31 percent of the treatment group completed at least one activity in either the fall or spring of the first year"
count if CompletedActivityFall16==1 | CompletedActivitySpring17 ==1
assert r(N)==1036
assert abs(r(N)/3300-.31)<0.005

* 740 completed the Spring 2017 wellness activities (out of 1059)
* "Over 75 percent of people who registered were signed up for Spring Into Motion"
count if CompletedActivitySpring17
assert r(N)==740

count if RegisteredActivitySpring17
assert r(N)==1059

count if ActivityS17=="Spring Into Motion"
assert r(N)/1059 >0.75

* 439 completed the Fall 2017 wellness activities (out of 771)
* "Over 70 percent of people who registered were signed up for Walktober"
count if CompletedActivityFall17
assert r(N)==439

count if RegisteredActivityFall17
assert r(N)==771

count if ActivityF17=="Walktober"
assert r(N)/771 >0.70


* 342 completed the Spring 2018 wellness activities (out of 607)
* "Over 60 percent of people who registered were signed up for Keep America Active"
count if CompletedActivitySpring18
assert r(N)==342

count if RegisteredActivitySpring18
assert r(N)==607

count if ActivityS18=="Keep America Active"
assert r(N)/607 >0.6

* "Twenty people who completed an activity were randomly selected to receive \$50 Amazon.com gift cards."
use "$Wellness_WhatDoesWWDo/data/raw/TreatmentAssignment-2017.dta", clear
count if Activity_50Dollar_Fall2017
assert r(N)==20
count if Activity_50Dollar_Spring2018
assert r(N)==20

**************************
* Claims data
**************************

* 8,323 people have claims between January 1, 2015 - July 31, 2017
* Appendix: "We obtained health insurance claims data for 8,461 university employees...who were members of Health Alliance at any point during the period January 1, 2015 through January 31, 2019."
* Note that "8,461" is also the number shown in Appendix Figure D.1
use "$Wellness_WhatDoesWWDo/data/raw/Claims-2019-03-31.dta", clear
keep if inrange(MonthOfService,ym(2015,01),ym(2019,01))
collapse (max) month_cov_frac, by(AnalysisID) fast
count if month_cov_frac>0
assert r(N) == 8461

* Appendix: "8,096 employees were members during the pre-period July 1, 2015 through July 31, 2016."
* Note: 8,096 should also match the sample size for Panel C in Balance Table 1b
use "$Wellness_WhatDoesWWDo/data/raw/Claims-2019-03-31.dta", clear
keep if inrange(MonthOfService,ym(2015,07),ym(2016,07))
collapse (max) month_cov_frac, by(AnalysisID) fast
count if month_cov_frac>0
assert r(N) == 8096

* Baseline spending histogram: "Data are from claims covering the period July 2015 - July 2016 ($N = 2,188$)"
use "$Wellness_WhatDoesWWDo/data/raw/Claims-2019-03-31.dta", clear
keep if inrange(MonthOfService,ym(2015,07),ym(2016,07))
merge m:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/raw/TreatmentAssignment-2016.dta", assert(match) nogenerate keepusing(StudyArm)
drop if inlist(StudyArm,"Not Randomized":StudyArm,"Control":StudyArm)
collapse (max) month_cov_frac, by(AnalysisID) fast
count if month_cov_frac>0
assert r(N) == 2188

* 339/3223= "11 percent of employees are not continuously enrolled throughout the 13-month pre-period".  294/3223= "9 percent are not continuously enrolled throughout the 12-month post-period"
use "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", clear
drop if mi(treat)
count if covg_0715_0716>0
assert r(N)==3223

count if covg_0715_0716>0 & covg_0715_0716<1
assert r(N)==339

count if covg_0816_0717>0 & covg_0816_0717<1
assert r(N)==294

***********************
* 2017 follow-up survey stats
***********************

* Out of the pre-registered set of 31 variables, 25 people are missing at least one response
local svy_hvars "everscreen active active_try cursmk othersmk formsmk drink drinkhvy chronic health1 health2 problems energy ehealth overweight badhealth sedentary"
local svy_uvars "druguse physician hospital"
local productivity_vars "sickdays hrsworked50 jobsatisf1 jobsatisf2 happywork presenteeism productive promotion jobsearch1 jobsearch2 mgmtsafety"
use "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", clear
keep if survey2017_c
gen miss=0
gen count = 0
qui foreach v in `svy_hvars' `svy_uvars' `productivity_vars' {
	replace miss=1 if mi(`v'_0717)
	replace count = count+1
}
count if miss==1
assert r(N)==25
assert count==31
drop miss count


* 4,824 people were sent invitations, 3,642 people clicked on survey, 3,611 began the survey, and 3,568 completed it
use "$Wellness_WhatDoesWWDo/data/raw/FollowupSurvey-2017.dta", clear
count if SentInvite
assert r(N)==4824

use "$Wellness_WhatDoesWWDo/data/raw/FollowupSurvey-2017.dta", clear
count if FinalDisposition_x10000=="Explicit refusal (withdrawn from study)":FinalDisposition_x10000
assert 4834-r(N)==4819

count if Clicked_FollowupSurvey==1
assert r(N)==3642

count if Completed_FollowupSurvey ==1
assert r(N)==3568
 
count if FinalDisposition_x10000=="Complete":FinalDisposition_x10000 | FinalDisposition_x10000=="Partial":FinalDisposition_x10000 | FinalDisposition_x10000=="Logged on to survey, not past consent language":FinalDisposition_x10000
assert r(N)==3642

count if FinalDisposition_x10000=="Complete":FinalDisposition_x10000 | FinalDisposition_x10000=="Partial":FinalDisposition_x10000
assert r(N)==3611

count if FinalDisposition_x10000=="Complete":FinalDisposition_x10000
assert r(N)==3568

* Completion rate for control: 75.4%, treat: 73.1% Treated group was 2.4% less likely to copmlete (p-val=0.079)
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", keep(match) assert(match using) keepusing(treat iThrive_talk_0717)
gen completed = FinalDisposition_x10000=="Complete":FinalDisposition_x10000
table treat, c(mean completed) format(%12.3fc)
reg completed treat, robust
assert abs(0.754-_b[_cons]) <0.001
assert abs(0.731-(_b[_cons]+_b[treat])) <0.001
assert abs(-0.024-_b[treat])<0.001
assert abs(tprob(e(df_r), abs(_b[treat]/_se[treat])) - 0.079) < 0.001

* Quickest completion time: 4.4 minutes
gen time_min = (TimeCompleted - TimeFirstClick )/(1000*60)
summ time_min

* Conditional on finishing in less than one hour, average completion time was 13.3 minutes
summ time_min if time_min<60
assert abs(13.3-r(mean))<0.1

* Out of 3,561 people who answered the age question, 20 were off by more than one year when compared to the HR data
keep if Completed_FollowupSurvey 
keep AnalysisID Background_Age treat iThrive_talk_0717 Time*
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/raw/HR_Baseline.dta", keep(match) assert(match using) nogenerate keepusing(Gender BirthYear)

count if !mi(Background_Age)
assert r(N)==3561

gen age_HR = 2017 - BirthYear
assert !mi(age_HR)
gen diff = age_HR - Background_Age

count if abs(diff)>1 & !mi(Background_Age)
assert r(N)==20

* 3% of control group, and 44% of treatment group, discussed iThrive with coworkers over past year.
table treat, c(mean iThrive_talk_0717) format(%12.2fc)
reg iThrive_talk_0717 ibn.treat, nocons
assert (_b[0.treat]-.0337079)<0.001
assert (_b[1.treat]-.4356846)<0.001

* "The August 2 reminder informed participants that ten people who completed the follow-up surveywould be chosen at random to receive a $100 Amazon.com gift card"
use "$Wellness_WhatDoesWWDo/data/raw/TreatmentAssignment-2017.dta", clear
count if Survey_100Dollar_2017
assert r(N)==10

***********************
* 2018 follow-up survey stats
***********************

* Out of the pre-registered set of 31 variables, 15 people are missing at least one response
local svy_hvars "everscreen active active_try cursmk othersmk formsmk drink drinkhvy chronic health1 health2 problems energy ehealth overweight badhealth sedentary"
local svy_uvars "druguse physician hospital"
local productivity_vars "sickdays hrsworked50 jobsatisf1 jobsatisf2 happywork presenteeism productive promotion jobsearch1 jobsearch2 mgmtsafety"
use "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", clear
keep if survey2018_c
gen miss=0
gen count = 0
qui foreach v in `svy_hvars' `svy_uvars' `productivity_vars' {
	replace miss=1 if mi(`v'_0718)
	replace count = count+1
}
count if miss==1
assert r(N)==15
assert count==31
drop miss count

* 4,800 people were sent email invitations, 3,120 people clicked on survey, and 3,020 completed it
use "$Wellness_WhatDoesWWDo/data/raw/FollowupSurvey-2018.dta", clear
count if SentInvite
assert r(N)==4800

count if Clicked_FollowupSurvey==1
assert r(N)==3120

count if Completed_FollowupSurvey ==1
assert r(N)==3020

* 15 people completed the 2018 survey after August 7
count if TimeCompleted > Cdhms(mdy(8,8,2018),1,1,1) & Completed_FollowupSurvey ==1
assert r(N)==15

* "Among those who completed the survey within an hour of clicking on the survey link for the first time, the average completion time was 16.9 minutes."
gen time_min = (TimeCompleted - TimeFirstClick )/(1000*60)
summ time_min if time_min<60
assert abs(16.9-r(mean))<.05

*"The completion rates for the control and treatment groups were 64.6 and 61.5 percent, respectively. This difference in completion rates is statistically significant ($p=0.036$)."
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", keep(match) assert(match using) keepusing(treat)
gen completed = FinalDisposition_x10000=="Complete":FinalDisposition_x10000
table treat, c(mean completed) format(%12.3fc)
reg completed treat, robust
assert abs(0.646-_b[_cons]) <0.001
assert abs(0.615-(_b[_cons]+_b[treat])) <0.001
assert abs(-0.031-_b[treat])<0.001
assert abs(tprob(e(df_r), abs(_b[treat]/_se[treat])) - 0.036) < 0.001

* Out of 3,009 people who answered the age question, 16 were off by more than one year when compared to the HR data
keep if Completed_FollowupSurvey 
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/raw/HR_Baseline.dta", keep(match) assert(match using) nogenerate keepusing(BirthYear)

count if !mi(Background_Age)
assert r(N)==3009

gen age_HR = 2018 - BirthYear
assert !mi(age_HR)
gen diff = age_HR - Background_Age

count if abs(diff)>1 & !mi(Background_Age)
assert r(N)==16

*****************************************************************
* Number of Participants at each stage by treatment/control group
*****************************************************************

use "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", clear

* Initial Pool = 12,459
count 
assert r(N) == 12459

* Study Sample = 4,834
count if StudyArm < 8
assert r(N) == 4834

* Control Group = 1,534
count if StudyArm == 7
assert r(N) == 1534

* Treatment Group = 3,300
count if StudyArm < 7
assert r(N) == 3300

* Treatment Biometric Screening + HRA Year 1 = 1,848
count if StudyArm < 7 & hra_c == 1
assert r(N) == 1848

* Treatment Fall Activities 2016 = 903
count if StudyArm < 7 & activity_f_c == 1
assert r(N) == 903

* Treatment Spring Activities 2016 = 740
count if StudyArm < 7 & activity_s_c == 1
assert r(N) == 740

* Control Follow Up Survey 2017 = 1,157
count if StudyArm == 7 & survey2017_c == 1
assert r(N) == 1157

* Treatment Follow Up Survey 2017 = 2,410
count if StudyArm < 7 & survey2017_c == 1
assert r(N) == 2410

* Control Screening 2017 = 595
count if StudyArm == 7 & screening2017_c == 1
assert r(N) == 595

* Treatment Screening + HRA 2017 = 1,272
count if StudyArm < 7 & hra_c_yr2 == 1
assert r(N) == 1272

* Treatment Fall Activities 2016 = 439
count if StudyArm < 7 & activity_f_c_yr2 == 1
assert r(N) == 439

* Treatment Spring Activities 2016 = 342
count if StudyArm < 7 & activity_s_c_yr2 == 1
assert r(N) == 342

* Control Survey 2018 = 991
count if StudyArm == 7 & survey2018_c == 1
assert r(N) == 991

* Treatment Survey 2018 = 2,029
count if StudyArm < 7 & survey2018_c == 1
assert r(N) == 2029

* Control Screening 2018 = 557
count if StudyArm == 7 & screening2018_c == 1
assert r(N) == 557

* Treatment Screening 2018 = 1,204
count if StudyArm < 7 & screening2018_c == 1
assert r(N) == 1204

* Treatment Groups A25 - C75 counts at each stage

local counts = "551 549 552 548 551 549"
local screen_hra = "246 270 305 339 335 353"
local well_fall = "102 167 115 197 133 189"
local well_spring = "82 132 94 174 99 159"

forval x = 1/6 {
	
	count if StudyArm == `x'
	assert r(N) == `:word `x' of `counts''
	
	count if StudyArm == `x' & screening_c == 1 & hra_c == 1
	assert r(N) == `:word `x' of `screen_hra''
	
	count if StudyArm == `x' & activity_f_c == 1
	assert r(N) == `:word `x' of `well_fall''
	
	count if StudyArm == `x' & activity_s_c == 1
	assert r(N) == `:word `x' of `well_spring''

}

* Follow-Up Sample = 4,834
count if StudyArm < 8
assert r(N) == 4834

* Online Survey 2017 = 3,567
count if survey2017_c == 1
assert r(N) == 3567

* Initial Pool for Year 2 == 4,819
count if StudyArm < 8 & StudyArm_2017 > 1
assert r(N) == 4819

* Online Survey 2017 = 3,567
count if survey2017_c == 1
assert r(N) == 3567

* Control Group = 1,530
count if StudyArm == 7 & StudyArm_2017 > 1
assert r(N) == 1530

* Treatment Group = 3,289
count if StudyArm < 7 & StudyArm_2017 > 1
assert r(N) == 3289

* Control groups C0 and C125
local counts = "765 765"
local screen = "192 403"

forval x = 1/2 {
	
	count if StudyArm_2017 == `x' + 3
	assert r(N) == `:word `x' of `counts''
	
	count if StudyArm_2017 == `x' + 3 & screening2017_c == 1
	assert r(N) == `:word `x' of `screen''

}


* Treatment Groups T0 and T125
local counts = "1644 1645"
local screen = "538 871"
local hra = "483 789"
local well_fall = "192 247"
local well_spring = "148 194"

forval x = 1/2 {
	
	count if StudyArm_2017 == `x' + 1
	assert r(N) == `:word `x' of `counts''
	
	count if StudyArm_2017 == `x' + 1 & screening2017_c == 1
	assert r(N) == `:word `x' of `screen''
	
	count if StudyArm_2017 == `x' + 1 & hra_c_yr2 == 1
	assert r(N) == `:word `x' of `hra''
	
	count if StudyArm_2017 == `x' + 1 & activity_f_c_yr2 == 1
	assert r(N) == `:word `x' of `well_fall''
	
	count if StudyArm_2017 == `x' + 1 & activity_s_c_yr2 == 1
	assert r(N) == `:word `x' of `well_spring''

}

* Second follow-up sample = 4,819
count if StudyArm < 8 & StudyArm_2017 > 1
assert r(N) == 4819

* Second follow-up survey = 3,020
count if survey2018_c == 1
assert r(N) == 3020

* Second screening = 1,761
count if screening2018_c == 1
assert r(N) == 1761

**************************************
* Number of Wellness Studies Ruled Out
**************************************

local lasso_dir = "$Wellness_WhatDoesWWDo/data/proc/0816_0717/lasso"
local winsorize_dir = "$Wellness_WhatDoesWWDo/results/intermediate_files/0816_0717/winsorization/estimates"
local estimate_dir = "$Wellness_WhatDoesWWDo/results/intermediate_files/0816_0717/treat_effects/estimates"


****
* ROI calcs
****

* "the lower bounds of the 99 percent confidence intervals for annual medical... are -\$396 ($=(17.2-2.577 \times 19.5)\times 12)$"
estimates use "`winsorize_dir'/ITT_Lasso_Winsorized_spendw01_0816_0717.ster"

assert (_b[rhs]-17.2)<0.05
assert (_se[rhs]-19.5)<0.05
assert abs((_b[rhs]-2.577*_se[rhs])*12 + 396)<0.5

* iThrive costs $152 per person, so we rule out ROI of 2.61 for spending
assert (396/152 - 2.61) < 0.005


* "we calculate that the lower bounds of the 99 percent confidence intervals for... absenteeism costs are... -\$91 ($=(0.138-2.577 \times 0.200) \times 240$)
estimates use "`estimate_dir'/ITT_Lasso_sickleave_0816_0717.ster"
assert (_b[rhs] -.138) <0.0005
assert (_se[rhs]-.200) <0.0005
assert abs((_b[rhs]-2.577*_se[rhs])*240 + 91)<0.5

* iThrive costs $152 per person, so we rule out ROI of .60 for absenteeism
assert (91/152 - .60) < 0.005

****
* Meta-analysis calcs
****
* Health ITT

estimates use "`winsorize_dir'/ITT_Lasso_Winsorized_spendw01_0816_0717.ster"

use "`lasso_dir'/lasso_data_spend_0816_0717.dta", clear
*sum spend_0816_0717 [aweight= covg_0816_0717] if in_lasso_IV == 1 & treat == 0

* Merge in winsorized spending
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", assert(match using) keep(match) nogen keepusing(spendw*)
assert _N==4834

sum spendw01_0816_0717 [aweight = covg_0816_0717] if in_lasso_IV == 1 & treat == 0

local health_itt = _b[rhs]/r(mean)
local health_itt_lb = (_b[rhs] - invttail(e(df_r),0.025)*_se[rhs])/r(mean)
local health_itt_ub = (_b[rhs] + invttail(e(df_r),0.025)*_se[rhs])/r(mean)

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "itt":estimate_type & outcome == "health":outcome
gen ruled_out = effect_percent < `health_itt_lb' | effect_percent > `health_itt_ub'
count
local n_health_itt = r(N)
count if ruled_out
local o_health_itt = r(N)

* Health IV

estimates use "`winsorize_dir'/IV_Lasso_Winsorized_spendw01_0816_0717.ster"

use "`lasso_dir'/lasso_data_spend_0816_0717.dta", clear
*sum spend_0816_0717 [aweight= covg_0816_0717] if in_lasso_IV == 1 & hra_c_nomiss == 1

* Merge in winsorized spending
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", assert(match using) keep(match) nogen keepusing(spendw*)
assert _N==4834

sum spendw01_0816_0717 [aweight= covg_0816_0717] if in_lasso_IV == 1 & hra_c_nomiss == 1

local health_iv = _b[rhs]/(r(mean) - _b[rhs])
local health_iv_lb = (_b[rhs] -  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])
local health_iv_ub = (_b[rhs] +  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "tot":estimate_type & outcome == "health":outcome
gen ruled_out = effect_percent < `health_iv_lb' | effect_percent > `health_iv_ub'
count
local n_health_iv = r(N)
count if ruled_out
local o_health_iv = r(N)


* Absenteeism ITT

estimates use "`estimate_dir'/ITT_Lasso_sickleave_0816_0717.ster"

use "`lasso_dir'/lasso_data_sickleave_0816_0717.dta", clear
sum sickleave_0816_0717 [aweight= sickdays_eligible_0816_0717] if in_lasso_IV == 1 & treat == 0

local absent_itt = _b[rhs]/r(mean)
local absent_itt_lb = (_b[rhs] - invttail(e(df_r),0.025)*_se[rhs])/r(mean)
local absent_itt_ub = (_b[rhs] + invttail(e(df_r),0.025)*_se[rhs])/r(mean)

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "itt":estimate_type & outcome == "absenteeism":outcome
gen ruled_out = effect_percent < `absent_itt_lb' | effect_percent > `absent_itt_ub'
count
local n_absent_itt = r(N)
count if ruled_out
local o_absent_itt = r(N)


* Absenteeism IV

estimates use "`estimate_dir'/IV_ztreat_Lasso_sickleave_0816_0717.ster"

use "`lasso_dir'/lasso_data_sickleave_0816_0717.dta", replace
sum sickleave_0816_0717 [aweight= sickdays_eligible_0816_0717] if in_lasso_IV == 1 & hra_c_nomiss == 1

local absent_iv = _b[rhs]/(r(mean) - _b[rhs])
local absent_iv_lb = (_b[rhs] -  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])
local absent_iv_ub = (_b[rhs] +  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "tot":estimate_type & outcome == "absenteeism":outcome
gen ruled_out = effect_percent < `absent_iv_lb' | effect_percent > `absent_iv_ub'
count
local n_absent_iv = r(N)
count if ruled_out
local o_absent_iv = r(N)


* "our 95-percent confidence interval rules out 20 of these 22 estimates."
assert `o_health_itt'==20
assert `n_health_itt'==22

* "our 95-percent confidence interval rules out 23 of the 33 studies."
assert `o_health_iv'==23
assert `n_health_iv'==33

* "Overall, our confidence intervals rule out 43 of 55 (78 percent) prior ITT and TOT point estimates for health spending."
assert abs( ((`o_health_itt'+`o_health_iv')/(`n_health_itt'+`n_health_iv')) - (43/55) ) < 0.001

* "our estimates rule out 51 of 57 (88 percent) prior ITT and TOT point estimates for absenteeism"
 assert abs( ((`o_absent_itt'+`o_absent_iv')/(`n_absent_itt'+`n_absent_iv')) - (51/57) ) < 0.001

* "Across both sets of outcomes, we rule out 94 of 112 (84 percent) prior estimates"
* Conclusion: "we can rule out 84 percent of the medical spending and absenteeism estimates from the prior literature"
* Intro: "Our 95 percent confidence intervals rule out 84 percent of the effects reported in 112 prior studies"
* Abstract: "Our 95\% confidence intervals rule out 84 percent of previous estimates on medical spending and absenteeism"
 assert abs( ((`o_health_itt'+`o_health_iv'+`o_absent_itt'+`o_absent_iv')/(`n_health_itt'+`n_health_iv'+`n_absent_itt'+`n_absent_iv')) - (94/112) ) < 0.001

di "Health ITT `o_health_itt'/`n_health_itt'"
di "Health IV `o_health_iv'/`n_health_iv'"
di "Total Health `=`o_health_itt'+`o_health_iv''/`=`n_health_itt'+`n_health_iv''"
di "Absent ITT `o_absent_itt'/`n_absent_itt'"
di "Absent IV `o_absent_iv'/`n_absent_iv'"
di "Total Absent `=`o_absent_itt'+`o_absent_iv''/`=`n_absent_itt'+`n_absent_iv''"
di "Total `=`o_health_itt'+`o_health_iv'+`o_absent_itt'+`o_absent_iv''/`=`n_health_itt'+`n_health_iv'+`n_absent_itt'+`n_absent_iv''"

******
* "If we do not winsorize medical spending, we rule out 40 of 55 (73 percent) prior health studies."
******
* Health ITT
estimates use "`estimate_dir'/ITT_Lasso_spend_0816_0717.ster"
use "`lasso_dir'/lasso_data_spend_0816_0717.dta", clear
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", assert(match using) keep(match) nogen keepusing(spendw*)
sum spend_0816_0717 [aweight = covg_0816_0717] if in_lasso_IV == 1 & treat == 0

local health_itt = _b[rhs]/r(mean)
local health_itt_lb = (_b[rhs] - invttail(e(df_r),0.025)*_se[rhs])/r(mean)
local health_itt_ub = (_b[rhs] + invttail(e(df_r),0.025)*_se[rhs])/r(mean)

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "itt":estimate_type & outcome == "health":outcome
gen ruled_out = effect_percent < `health_itt_lb' | effect_percent > `health_itt_ub'
count if ruled_out
local o_health_itt = r(N)

* Health IV
estimates use "`estimate_dir'/IV_ztreat_Lasso_spend_0816_0717.ster"
use "`lasso_dir'/lasso_data_spend_0816_0717.dta", clear
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", assert(match using) keep(match) nogen keepusing(spendw*)
assert _N==4834
sum spend_0816_0717 [aweight= covg_0816_0717] if in_lasso_IV == 1 & hra_c_nomiss == 1

local health_iv = _b[rhs]/(r(mean) - _b[rhs])
local health_iv_lb = (_b[rhs] -  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])
local health_iv_ub = (_b[rhs] +  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "tot":estimate_type & outcome == "health":outcome
gen ruled_out = effect_percent < `health_iv_lb' | effect_percent > `health_iv_ub'
count if ruled_out
local o_health_iv = r(N)

assert `o_health_itt'+`o_health_iv'==40
assert `n_health_itt'+`n_health_iv'==55

di "Total Health `=`o_health_itt'+`o_health_iv''/`=`n_health_itt'+`n_health_iv''"


******
* "if we restrict our comparison to only the set of RCTs, we rule out 21 of 22 
* (95 percent) prior estimates."
******

* Health ITT

estimates use "`winsorize_dir'/ITT_Lasso_Winsorized_spendw01_0816_0717.ster"

use "`lasso_dir'/lasso_data_spend_0816_0717.dta", clear
*sum spend_0816_0717 [aweight= covg_0816_0717] if in_lasso_IV == 1 & treat == 0

* Merge in winsorized spending
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", assert(match using) keep(match) nogen keepusing(spendw*)
assert _N==4834

sum spendw01_0816_0717 [aweight = covg_0816_0717] if in_lasso_IV == 1 & treat == 0

local health_itt = _b[rhs]/r(mean)
local health_itt_lb = (_b[rhs] - invttail(e(df_r),0.025)*_se[rhs])/r(mean)
local health_itt_ub = (_b[rhs] + invttail(e(df_r),0.025)*_se[rhs])/r(mean)

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "itt":estimate_type & outcome == "health":outcome & study_type == "RCT":study_type
gen ruled_out = effect_percent < `health_itt_lb' | effect_percent > `health_itt_ub'
count
local n_health_itt = r(N)
count if ruled_out
local o_health_itt = r(N)

* Health IV

estimates use "`winsorize_dir'/IV_Lasso_Winsorized_spendw01_0816_0717.ster"

use "`lasso_dir'/lasso_data_spend_0816_0717.dta", clear
*sum spend_0816_0717 [aweight= covg_0816_0717] if in_lasso_IV == 1 & hra_c_nomiss == 1

* Merge in winsorized spending
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", assert(match using) keep(match) nogen keepusing(spendw*)
assert _N==4834

sum spendw01_0816_0717 [aweight= covg_0816_0717] if in_lasso_IV == 1 & hra_c_nomiss == 1

local health_iv = _b[rhs]/(r(mean) - _b[rhs])
local health_iv_lb = (_b[rhs] -  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])
local health_iv_ub = (_b[rhs] +  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "tot":estimate_type & outcome == "health":outcome & study_type == "RCT":study_type
gen ruled_out = effect_percent < `health_iv_lb' | effect_percent > `health_iv_ub'
count
local n_health_iv = r(N)
count if ruled_out
local o_health_iv = r(N)

* Absenteeism ITT

estimates use "`estimate_dir'/ITT_Lasso_sickleave_0816_0717.ster"

use "`lasso_dir'/lasso_data_sickleave_0816_0717.dta", clear
sum sickleave_0816_0717 [aweight= sickdays_eligible_0816_0717] if in_lasso_IV == 1 & treat == 0

local absent_itt = _b[rhs]/r(mean)
local absent_itt_lb = (_b[rhs] - invttail(e(df_r),0.025)*_se[rhs])/r(mean)
local absent_itt_ub = (_b[rhs] + invttail(e(df_r),0.025)*_se[rhs])/r(mean)

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "itt":estimate_type & outcome == "absenteeism":outcome & study_type == "RCT":study_type
gen ruled_out = effect_percent < `absent_itt_lb' | effect_percent > `absent_itt_ub'
count
local n_absent_itt = r(N)
count if ruled_out
local o_absent_itt = r(N)


* Absenteeism IV

estimates use "`estimate_dir'/IV_ztreat_Lasso_sickleave_0816_0717.ster"

use "`lasso_dir'/lasso_data_sickleave_0816_0717.dta", replace
sum sickleave_0816_0717 [aweight= sickdays_eligible_0816_0717] if in_lasso_IV == 1 & hra_c_nomiss == 1

local absent_iv = _b[rhs]/(r(mean) - _b[rhs])
local absent_iv_lb = (_b[rhs] -  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])
local absent_iv_ub = (_b[rhs] +  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "tot":estimate_type & outcome == "absenteeism":outcome & study_type == "RCT":study_type
gen ruled_out = effect_percent < `absent_iv_lb' | effect_percent > `absent_iv_ub'
count
local n_absent_iv = r(N)
count if ruled_out
local o_absent_iv = r(N)

assert abs( ((`o_health_itt'+`o_health_iv'+`o_absent_itt'+`o_absent_iv')/(`n_health_itt'+`n_health_iv'+`n_absent_itt'+`n_absent_iv')) - (21/22) ) < 0.001

di "Health ITT `o_health_itt'/`n_health_itt'"
di "Health IV `o_health_iv'/`n_health_iv'"
di "Total Health `=`o_health_itt'+`o_health_iv''/`=`n_health_itt'+`n_health_iv''"
di "Absent ITT `o_absent_itt'/`n_absent_itt'"
di "Absent IV `o_absent_iv'/`n_absent_iv'"
di "Total Absent `=`o_absent_itt'+`o_absent_iv''/`=`n_absent_itt'+`n_absent_iv''"
di "Total `=`o_health_itt'+`o_health_iv'+`o_absent_itt'+`o_absent_iv''/`=`n_health_itt'+`n_health_iv'+`n_absent_itt'+`n_absent_iv''"

******
* "If we restrict our comparison to just the studies that lasted 12 months or 
* less, we rule out 39 of 47 (83 percent) prior estimates"
******

* Health ITT

estimates use "`winsorize_dir'/ITT_Lasso_Winsorized_spendw01_0816_0717.ster"

use "`lasso_dir'/lasso_data_spend_0816_0717.dta", clear
*sum spend_0816_0717 [aweight= covg_0816_0717] if in_lasso_IV == 1 & treat == 0

* Merge in winsorized spending
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", assert(match using) keep(match) nogen keepusing(spendw*)
assert _N==4834

sum spendw01_0816_0717 [aweight = covg_0816_0717] if in_lasso_IV == 1 & treat == 0

local health_itt = _b[rhs]/r(mean)
local health_itt_lb = (_b[rhs] - invttail(e(df_r),0.025)*_se[rhs])/r(mean)
local health_itt_ub = (_b[rhs] + invttail(e(df_r),0.025)*_se[rhs])/r(mean)

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "itt":estimate_type & outcome == "health":outcome & timelength <= 12
gen ruled_out = effect_percent < `health_itt_lb' | effect_percent > `health_itt_ub'
count
local n_health_itt = r(N)
count if ruled_out
local o_health_itt = r(N)

* Health IV

estimates use "`winsorize_dir'/IV_Lasso_Winsorized_spendw01_0816_0717.ster"

use "`lasso_dir'/lasso_data_spend_0816_0717.dta", clear
*sum spend_0816_0717 [aweight= covg_0816_0717] if in_lasso_IV == 1 & hra_c_nomiss == 1

* Merge in winsorized spending
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", assert(match using) keep(match) nogen keepusing(spendw*)
assert _N==4834

sum spendw01_0816_0717 [aweight= covg_0816_0717] if in_lasso_IV == 1 & hra_c_nomiss == 1

local health_iv = _b[rhs]/(r(mean) - _b[rhs])
local health_iv_lb = (_b[rhs] -  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])
local health_iv_ub = (_b[rhs] +  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "tot":estimate_type & outcome == "health":outcome & timelength <= 12
gen ruled_out = effect_percent < `health_iv_lb' | effect_percent > `health_iv_ub'
count
local n_health_iv = r(N)
count if ruled_out
local o_health_iv = r(N)


* Absenteeism ITT

estimates use "`estimate_dir'/ITT_Lasso_sickleave_0816_0717.ster"

use "`lasso_dir'/lasso_data_sickleave_0816_0717.dta", clear
sum sickleave_0816_0717 [aweight= sickdays_eligible_0816_0717] if in_lasso_IV == 1 & treat == 0

local absent_itt = _b[rhs]/r(mean)
local absent_itt_lb = (_b[rhs] - invttail(e(df_r),0.025)*_se[rhs])/r(mean)
local absent_itt_ub = (_b[rhs] + invttail(e(df_r),0.025)*_se[rhs])/r(mean)

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "itt":estimate_type & outcome == "absenteeism":outcome & timelength <= 12
gen ruled_out = effect_percent < `absent_itt_lb' | effect_percent > `absent_itt_ub'
count
local n_absent_itt = r(N)
count if ruled_out
local o_absent_itt = r(N)

* Absenteeism IV

estimates use "`estimate_dir'/IV_ztreat_Lasso_sickleave_0816_0717.ster"

use "`lasso_dir'/lasso_data_sickleave_0816_0717.dta", replace
sum sickleave_0816_0717 [aweight= sickdays_eligible_0816_0717] if in_lasso_IV == 1 & hra_c_nomiss == 1

local absent_iv = _b[rhs]/(r(mean) - _b[rhs])
local absent_iv_lb = (_b[rhs] -  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])
local absent_iv_ub = (_b[rhs] +  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "tot":estimate_type & outcome == "absenteeism":outcome & timelength <= 12
gen ruled_out = effect_percent < `absent_iv_lb' | effect_percent > `absent_iv_ub'
count
local n_absent_iv = r(N)
count if ruled_out
local o_absent_iv = r(N)


assert abs( ((`o_health_itt'+`o_health_iv'+`o_absent_itt'+`o_absent_iv')/(`n_health_itt'+`n_health_iv'+`n_absent_itt'+`n_absent_iv')) - (39/47) ) < 0.001


di "Health ITT `o_health_itt'/`n_health_itt'"
di "Health IV `o_health_iv'/`n_health_iv'"
di "Total Health `=`o_health_itt'+`o_health_iv''/`=`n_health_itt'+`n_health_iv''"
di "Absent ITT `o_absent_itt'/`n_absent_itt'"
di "Absent IV `o_absent_iv'/`n_absent_iv'"
di "Total Absent `=`o_absent_itt'+`o_absent_iv''/`=`n_absent_itt'+`n_absent_iv''"
di "Total `=`o_health_itt'+`o_health_iv'+`o_absent_itt'+`o_absent_iv''/`=`n_health_itt'+`n_health_iv'+`n_absent_itt'+`n_absent_iv''"

******
* "If wecombine RCTs and studies that use a pre/post design, 
* we continue to rule out 68 of 81 (84 percent) prior estimates."
******

* Health ITT

estimates use "`winsorize_dir'/ITT_Lasso_Winsorized_spendw01_0816_0717.ster"

use "`lasso_dir'/lasso_data_spend_0816_0717.dta", clear
*sum spend_0816_0717 [aweight= covg_0816_0717] if in_lasso_IV == 1 & treat == 0

* Merge in winsorized spending
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", assert(match using) keep(match) nogen keepusing(spendw*)
assert _N==4834

sum spendw01_0816_0717 [aweight = covg_0816_0717] if in_lasso_IV == 1 & treat == 0

local health_itt = _b[rhs]/r(mean)
local health_itt_lb = (_b[rhs] - invttail(e(df_r),0.025)*_se[rhs])/r(mean)
local health_itt_ub = (_b[rhs] + invttail(e(df_r),0.025)*_se[rhs])/r(mean)

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "itt":estimate_type & outcome == "health":outcome & (study_type == "RCT":study_type | pre_post)
gen ruled_out = effect_percent < `health_itt_lb' | effect_percent > `health_itt_ub'
count
local n_health_itt = r(N)
count if ruled_out
local o_health_itt = r(N)

* Health IV

estimates use "`winsorize_dir'/IV_Lasso_Winsorized_spendw01_0816_0717.ster"

use "`lasso_dir'/lasso_data_spend_0816_0717.dta", clear
*sum spend_0816_0717 [aweight= covg_0816_0717] if in_lasso_IV == 1 & hra_c_nomiss == 1

* Merge in winsorized spending
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", assert(match using) keep(match) nogen keepusing(spendw*)
assert _N==4834

sum spendw01_0816_0717 [aweight= covg_0816_0717] if in_lasso_IV == 1 & hra_c_nomiss == 1

local health_iv = _b[rhs]/(r(mean) - _b[rhs])
local health_iv_lb = (_b[rhs] -  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])
local health_iv_ub = (_b[rhs] +  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "tot":estimate_type & outcome == "health":outcome & (study_type == "RCT":study_type | pre_post)
gen ruled_out = effect_percent < `health_iv_lb' | effect_percent > `health_iv_ub'
count
local n_health_iv = r(N)
count if ruled_out
local o_health_iv = r(N)


* Absenteeism ITT

estimates use "`estimate_dir'/ITT_Lasso_sickleave_0816_0717.ster"

use "`lasso_dir'/lasso_data_sickleave_0816_0717.dta", clear
sum sickleave_0816_0717 [aweight= sickdays_eligible_0816_0717] if in_lasso_IV == 1 & treat == 0

local absent_itt = _b[rhs]/r(mean)
local absent_itt_lb = (_b[rhs] - invttail(e(df_r),0.025)*_se[rhs])/r(mean)
local absent_itt_ub = (_b[rhs] + invttail(e(df_r),0.025)*_se[rhs])/r(mean)

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "itt":estimate_type & outcome == "absenteeism":outcome & (study_type == "RCT":study_type | pre_post)
gen ruled_out = effect_percent < `absent_itt_lb' | effect_percent > `absent_itt_ub'
count
local n_absent_itt = r(N)
count if ruled_out
local o_absent_itt = r(N)

* Absenteeism IV

estimates use "`estimate_dir'/IV_ztreat_Lasso_sickleave_0816_0717.ster"

use "`lasso_dir'/lasso_data_sickleave_0816_0717.dta", replace
sum sickleave_0816_0717 [aweight= sickdays_eligible_0816_0717] if in_lasso_IV == 1 & hra_c_nomiss == 1

local absent_iv = _b[rhs]/(r(mean) - _b[rhs])
local absent_iv_lb = (_b[rhs] -  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])
local absent_iv_ub = (_b[rhs] +  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "tot":estimate_type & outcome == "absenteeism":outcome & (study_type == "RCT":study_type | pre_post)
gen ruled_out = effect_percent < `absent_iv_lb' | effect_percent > `absent_iv_ub'
count
local n_absent_iv = r(N)
count if ruled_out
local o_absent_iv = r(N)

assert abs( ((`o_health_itt'+`o_health_iv'+`o_absent_itt'+`o_absent_iv')/(`n_health_itt'+`n_health_iv'+`n_absent_itt'+`n_absent_iv')) - (68/81) ) < 0.001


di "Health ITT `o_health_itt'/`n_health_itt'"
di "Health IV `o_health_iv'/`n_health_iv'"
di "Total Health `=`o_health_itt'+`o_health_iv''/`=`n_health_itt'+`n_health_iv''"
di "Absent ITT `o_absent_itt'/`n_absent_itt'"
di "Absent IV `o_absent_iv'/`n_absent_iv'"
di "Total Absent `=`o_absent_itt'+`o_absent_iv''/`=`n_absent_itt'+`n_absent_iv''"
di "Total `=`o_health_itt'+`o_health_iv'+`o_absent_itt'+`o_absent_iv''/`=`n_health_itt'+`n_health_iv'+`n_absent_itt'+`n_absent_iv''"

*********************************
* Publication Bias Calculations *
*********************************

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear

keep if has_se
drop if drop_dup

keep pub_effect pub_se
export delimited using "$Wellness_WhatDoesWWDo/data/proc/lit_review/pub_bias.csv", novarnames replace

preserve
* .csv file is then uploaded to online interface at https://maxkasy.github.io/home/metastudy/
* Results for publication analysis is return on website, using a symmetric distribution for
* p, a cutoff for p at 1.96, and a Student-t distribution for the model
* The following lines of code take those results and makes a .tex table

* theta, tau, nu, and beta point estimates (produced at the website above):
matrix define pub_bias = (-0.583, 0.385, 1.990, 0.369\0.398, 0.302, 0.478, 0.153)

* Store and format these estimates
drop _all
svmat pub_bias
tostring *, format(%12.3fc) force replace
forval x = 1/4 {
	replace pub_bias`x' = "(" + pub_bias`x' + ")" if _n==2
}

* Create table
local headerlines `"{(1)}&{(2)}&{(3)}&{(4)}" "\midrule & & & Relative Publication" "\multicolumn{3}{c}{Student-\(t\) Distribution of Wellness Program Effects} & Probability for \(p>0.05\)" "\cmidrule(lr){1-3} \cmidrule(lr){4-4} \(\bar{\theta}\) & \(\tilde{\tau}\) & \(\nu\) & \(\beta_{1,p}\)"'
local fn "Notes: Using the method of \cite{Andrews:2019}, we estimate the model given by Equation \ref{eqn:pub_bias}. Table reports the meta-analysis estimates of the bias-corrected distribution of wellness program effects on medical spending and absenteeism, with location (\(\bar{\theta}\)), scale (\(\tilde{\tau}\)), and degrees of freedom (\(\tilde{\nu}\)) parameters for a Student-\(t\) distribution. Publication probability \(\beta_{1,p}\) for studies not significant at the 5\% level is normalized relative to studies significant at the 5\% level."
texsave using "$Wellness_WhatDoesWWDo/results/tables/appendix_pub_bias.tex", headerlines("`headerlines'") title("Meta-Analysis Estimates of Publication Bias") footnote(`"`fn'"', size(footnotesize)) marker(tab:appendix_pub_bias) replace nonames nofix
restore

***
* "The bias-corrected mean estimate is not significantly difference from zero (p=0.14)."
***

assert abs(2*(1-normal(abs(-0.583/0.398))) - 0.14) < 0.05

***
* "For ease of presentation, we omit 3 studies from this figure with extreme values."
***

count
local count0 = r(N)

local trim = 1000

drop if abs(pub_effect) > `trim'
drop if pub_se > `trim'/1.96 & pub_se != .

count
local count1 = r(N)

assert abs(`count1'/`count0' - 37/40) < 0.001

** EOF
